DUDLEY  KNOX  LIBRARY 
MAVAL  POSTGRADUATE  SCHOOL 
MONTEREY    CALIFORNIA  93943-500', 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


THESIS 

AN  IMPLEMENTATION 

OF  THE 

PROJECTIVE  ALGORITHM  FOR  LINEAR  PROGRAMMING 

by 

Guenter  W.  Bretschneider 

September  1985 

Thesis 

Advisor:                  R. 

Kevin 

Wood 

Approved  for  public  release;  distribution  is  unlimited 


T226039 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whan  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


1.  REPORT  NUMBER 


2.  GOVT  ACCESSION  NO 


3.   RECIPIENT'S  CATALOG  NUMBER 


4.     TITLE  (and  Subtitle) 

An  Implementation  of  the  Projective 
Algorithm  for  Linear  Programming 


5.     TYPE  OF   REPORT  a   PERIOD  COVERED 

Master's  Thesis 
September,  1985 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHORCs; 

Guenter  W.  Bretschneider 


8.  CONTRACT  OR  GRANT  NUMBER(«; 


9.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California   93943-5100 


10.     PROGRAM  ELEMENT.  PROJECT,   TASK 
AREA  4   WORK   UNIT  NUMBERS 


II.     CONTROLLING  OFFICE   NAME   AND   ADDRESS 

Naval  Postgraduate  School 
Monterey,  California   93943-5100 


12.     REPORT   DATE 

September,    19  8  5 


13.     NUMBER  OF   PAGES 

43 


U.     MONITORING  AGENCY  NAME  4    AODRESSf//  different  from  Controlling  Office) 


15.     SECURITY   CLASS,  (of  thla  report) 


15«.     DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


16.     DISTRIBUTION   STATEMENT  rot  this  Report) 


Approved  for  public  release;  distribution  is  unlimited 


17.     DISTRIBUTION  STATEMENT  (oi  the  abstract  entered  In  Block  30,   If  different  from  Report) 


18.     SUPPLEMENTARY   NOTES 


19.     KEY   WORDS  r Continue  on  reverse   side  it  necessary  and  Identity  by  block  number) 


Linear  Programming,  Projection  Method,  Symmetric  Matrices, 
Least  Squares  Problem,  Cholesky  Factorization 


20.      ABSTRACT  f Continue  on  reverse   aide   If  naceaaary  and  Identity  by  block  number) 

An  algorithm  to  solve  linear  programming  problems  is  presented 
which  is  based  on  Karmarkar's  projective  method.   The  algorithm 
includes  a  practical  method  to  project  a  general  linear  program- 
ming problem  onto  a  unit  simplex  and  eliminates  the  a  priori  need 
to  know  the  optimal  value  of  the  objective  function.   The  imple- 
mentation conserves  sparsity.   The  key  part  of  the  implementation 
is  the  solution  of  a  linear  least-squares  problem  to  find  an 


DD     |    JAN    73     1473  EDITION   OF    1   NOV  65  IS  OBSOLETE 

S    N   0102-  LF-  014-  6601  1 


SECURITY   CLASSIFICATION   OF   THIS  PAGE  (When  Data  Bntaraa) 


SECURITY   CLASSIFICATION  OF   THIS  PAGE  (Whmn  Dmm  Bnft»d) 


20.   improving  direction:   a  direct  and  an  iterative  method  are 

implemented  to  solve  this  problem.   The  direct  method  employs 
the  minimum-degree  heuristic  to  reorder  the  system  of  normal 
equations,  and  thus  conserve  sparsity  during  the  following 
Cholesky  factor  of  the  normal  equation  martrix  as  a  precon- 
ditioner  for  conjugate  gradient  iterations  which  are  per- 
formed implicitly  on  the  preconditioned  matrix.   The  study 
concludes  with  implementation  remarks,  and  computational 
results . 


SECURITY   CLASSIFICATION   OF   THIS  P  AGZ(Wh»n  Dmlm  Ent;,d) 


Approved  for  public  release;    distribution  is  unlimited 


An  Implementation 

of  the 

Projective  Algorithm  for  Linear  Programming" 


by 


Guenter  W.    Bretschneider 

Captain,    German  Air  Force 

M.S.    in  Aeronautical  Engineering, 

Hochschule  der  Bundeswehr,    Munich,    1976 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
September   1985 


ABSTRACT 

An  algorithm  to  solve  linear  programming  problems  is  presented 
which  is  based  on  Karmarkar's  projective  method.  The  algorithm 
includes  a  practical  method  to  project  a  general  linear  programming 
problem  onto  a  unit  simplex  and  eliminates  the  a  priori  need  to  know 
the  optimal  value  of  the  objective  function.  The  implementation 
conserves  sparsity.  The  key  part  of  the  implementation  is  the  solution 
of  a  linear  least- squares  problem  to  find  an  improving  direction:  a 
direct  and  an  iterative  method  are  implemented  to  solve  this  problem. 
The  direct  method  employs  the  minimum- degree  heuristic  to  reorder  the 
system  of  normal  equations ,  and  thus  conserve  sparsity  during  the 
following  Cholesky  factorization.  The  iterative  method  uses  the  incom- 
plete Cholesky  factor  of  the  normal  equation  matrix  as  a  preconditioner 
for  conjugate  gradient  iterations  which  are  performed  implicitly  on  the 
preconditioned  matrix.  The  study  concludes  with  implementation 
remarks,    and  computational  results. 


Dudley  kno:       ^haky 

NAVAL  POSTGitis-jJliATB  SCHOOL 
MONTEREY,  CALIFORNIA  93943-5002 


TABLE  OF  CONTENTS 


I.  INTRODUCTION 8 

A.  THE  BASIC  PROJECTIVE  METHOD 8 

B.  DERIVATION  OF  THE  ALGORITHM       10 

1 .  The   Canonical  Form       10 

2.  The  Projective  Transformation       11 

3.  Optimization  Over  a  Sphere 12 

II.  •         A  VARIANT  OF  THE  PROJECTIVE  METHOD       17 

A.  INTRODUCTION 17 

B.  DERIVATION  OF  THE  MODIFIED  ALGORITHM 17 

C.  THE  LINEAR  LEAST- SQUARES  PROBLEM       19 

D.  IMPLEMENTATION 21 

III.  SOLVING  THE  LINEAR  LEAST- SQUARES  PROBLEM 26 

A.  INTRODUCTION 26 

B.  CHOLESKY  FACTORIZATION       26 

1.  The  Minimum- Degree   Ordering  Heuristic        27 

2 .  Factorization  and   Solution        28 

C.  INCOMPLETE   CHOLESKY  FACTORIZATION 29 

1.  The   Incomplete   Factorization 31 

2 .  The   Conjugate  Gradient  Iteration 33 

IV.  MODIFICATIONS   AND   COMPUTATIONAL  RESULTS      35 

A.  ALGORITHMIC  MODIFICATIONS 35 

1.  Iterative   Improvements 35 

2 .  Removal  of  the  Artificial  Column       35 

T 

3.  Positive   Semidefinite   BB       36 

4.  Weighted  Homogeneity  Variable 36 

B.  TEST  PROBLEMS   AND   COMPUTATIONAL  RESULTS        .    .  37 

C.  CONCLUSIONS       40 

LIST  OF  REFERENCES 41 


INITIAL  DISTRIBUTION  LIST       43 


LIST  OF  TABLES 

1.  ALGORITHM  1 16 

2.  ALGORITHM  2 20 

3.  MINIMUM-DEGREE  ORDERING  /  CHOLESKY 
FACTORIZATION  ALGORITHM   (MDOC) 27 

4.  INCOMPLETE  CHOLESKY  -   CONJUGATE  GRADIENT 

(ICCG)   ALGORITHM 31 

5.  TEST  PROBLEMS 37 

6.  •         COMPUTATIONAL  RESULTS       38 

7.  DENSITIES 39 


I.    INTRODUCTION 

In  recent  papers  Karmarkar  [Refs.  1,2]  has  presented  a  new 
method  for  the  solution  of  linear  programming  (LP)  problems;  His  new 
solution  technique  moves  from  some  feasible  starting  point  across  the 
interior  region  of  a  polytope  that  is  defined  by  the  problem  constraints. 
He  shows  that  the  number  of  steps  to  find  an  optimal  solution  with  his 
technique  is  polynomially  bounded. 

In  contrast,  the  simplex  algorithm  which,  is  widely  used  for  the 
solution  of  LP  problems  finds  an  optimal  solution  by  moving  from  vertex 
to  vertex  on  the  polytope.  This  is  known,  in  the  worst  case,  to  require 
an  exponential  number  of  steps    [Ref .    3] . 

The  polynomial  bound  of  the  projective  algorithm  makes  the  new 
solution  technique  very  appealing  to  researchers.  The  theory  of  the 
new  algorithm  seems  to  be  widely  accepted  among  experts,  while 
Karmarkar 's  claim  that  his  algorithm  is  50-100  times  faster  than  the 
simplex  method  has  met   skepticism   [Ref.    4] . 

In  this  study  a  variant  of  the  projective  method  is  implemented, 
and   some  well  known  test  problems  are   solved. 

A.       THE  BASIC  PROJECTIVE  METHOD 

Following  Karmarkar  [Ref.  1:  p.  4],  the  number  of  steps  of  the 
algorithm  depends  on  R/r,  where  R  is  the  radius  of  the  sphere 
circumscribing  the  polytope,    and  r  the  radius  of  the  inscribed   sphere. 

Assume  a   general  LP  of  the  form 


Min 

T 

s.t. 

A  x 

= 

b 

X 

> 

0 

(LP1) 


With   the   assumption   that   the   sum   of   its   variables   has   an   upper   bound, 
and  with  the  proper  scaling  of  variables,    a  convexity  constraint 


lTx  =  1  (l.l) 


can  be  added  to  LP1.  This  transformation  maps  the  LP  onto  a  unit 
simplex  S  whose  center  is  at  a  =  (1/n,  .  .  .  ,  1/n) .  Then,  a  transforma- 
tion is  performed  such  that  a  feasible  interior  starting  point  is  mapped 
onto  a   . 

If  B(a  , r)  is  the  largest  sphere  with  center  a  that  can  be 
inscribed  into  the  simplex  S,  and  B(a  , R)  is  the  smallest  sphere 
circumscribing   S,    then 

R/r  =   n.  (1.2) 

By  restricting  a  solution  x  to  remain  inside  the  largest  inscribed  sphere 
B(a  ,  r) ,  the  method  achieves  in  one  iteration  a  reduction  in  the  differ- 
ence between  the  current  objective  value  and  the  optimal  objective  value 
by  a  factor  of  (1  -  1/n).  Following  Shanno  [Ref.  5]  a  simple  proof  is: 
Let 

f*  =  min  cTx,        x  6   S,    Ax  -  b,  (1.3) 

f     =  min  cTx,        x  €  B(aQ,r),    Ax  =  b,  lTx  =   1,                                (1.4) 

f     =  min  cTx,        x  e  B(aQ,R),    Ax  =  b,  lTx  =   1.                               (1.5) 

Then, 

cTa0  -  1    <     cTa0  -  r     <     cTa0  -  f ,  (1.6) 

and  with  equation   1.2 

cTa0  -   f     =     n(cTa0  -  _f )  .  (1.7) 

From  that 

(f  -   r)/(cTa0  -   T)   <    (1   -    1/n).  (1.8) 


If  a  linear  objective  function  could  be  maintained  at  each  iteration, 
it  follows  that  an  upper  bound  on  the  number  of  steps  required  to  find 
an  optimal  solution  is  of  O(n).  Unfortunately,  the  projective  transfor- 
mations needed  to  continue  the  algorithm  result  in  a  nonlinear  objective 
function. 

B.       DERIVATION  OF  THE  ALGORITHM 

1.       The  Canonical  Form 

Suppose  we  have  a  linear  programming  problem  of  the  form 
LP1.  Karmarkar's  method  requires  that  this  LP  be  transformed  into  the 
following  canonical  form, 


Min 

T 
ex 

s.t. 

A  x 

= 

0 

lTx 

= 

1 

X 

> 

0 

(LP2) 


where  the  optimal  solution  value  is  zero.  With  the  assumption  that  the 
sum  of  the  variables  of  an  LP  is  bounded  above  and  subsequent  scaling 
by  this  bound,  a  convexity  constraint  (equation  1.1)  can  be  added  to 
LP1.  In  practice  this  can  be  a  problem  because  choosing  too  large  an 
upper  bound  may  cause  numerical  problems. 

LP2     requires     that     the     nonhomogeneous     system    of     equations 

Ax  :  b    be    made    homogeneous.     Karmarkar     [Ref.    1:     p. 34]     proposes    a 

T 
transformation  that  would  transform  the  i-th  equation  a-    x  =   b-      to 

2(ajj   -   bi)x]  =  0.  (1.9) 

This  transformation  has  the  disadvantage  that  when  b  is  dense  the 
sparsity  of  A  will  be  lost.  Karmarkar  requires  the  optimal  solution 
value  to  LP2  to  be  zero  together  with  a  special  stopping  criterion 
(equation  1.26),  to  prove  the  polynomial  bound.  To  achieve  the  zero 
objective      value      the      optimal      objective      function      value      f        of      the 
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untransformed  LP  has   to  be  known  in  advance,    and   the  following   trans- 
formation    made 

cTx  -   £*  =  cTx  -  f*lTx  =    (c  -ri)Tx   .  (1.10) 


Karmarkar  concludes  [Ref.  2:  p.  387]  that  if  the  minimal  objective  value 
determined  by  the  projective  algorithm  is  not  equal  to  zero,  the  original 
problem  must  be  either  infeasible  or  unbounded.  Another  method  for 
making  the  transformation  from  LP1  to  LP2  is  discussed  in 
Chapter  II.  B. 

2.       The  Projective  Transformation 

Let  x     >   0  be  a  known  feasible   solution   to  LP2 .    The   invertible 
projective  transformation 


y  = 


D_1x  (1.11) 


lTD_1x 


T 

maps     any     x,     such    that       1    x     =     1       and       x  >  0,     onto     y     such    .that 

T 

1    y  =   1      and     y  >  0.        D     is     a      diagonal     matrix     whose     entries     are 

(x2    >  ■  ■  ■  >  xn   )  • 

The   transformation   maps    the   LP2    unit   simplex   in   x-  space   onto 

another     unit      simplex      in     y-  space.      The     point     x        is     mapped     onto 

y°   =    1/n   1     ,    the   center  of   the   unit    simplex   in   y-space.    The    inverse   of 

the   transformation    (equation   1.11)    is   given  by 


D  y  (1.12) 

x  =   . 

1TD  y 


After  the  transformation,    LP2   can  now  be  restated  as 
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cTD  y 

Min     


1TD  y 


s.t.      A  D  y      =0  (LP3) 

=     1 
>     0   . 


iT  y     =    i 


Strictly  speaking  LP3  is  not  a  linear  program  because  its  objective 
function  is  a  rational  function  of  y.  Karmarkar  [Ref.  1:  page  18ff] 
shows  that  it  is  sufficient  to  consider  a  linear  approximation  to  the 
objective  function  of  LP3,    which  leads  to 


Min 

cxD  y 

s.t. 

A  D  y 

=     0 

iTy 

=    i 

y 

>     0 

(LP4) 


3.       Optimization  Over  a  Sphere 

A   solution   to   LP4   is   now   restricted   to  lie   within   a   sphere   with 
center  at  y°  and  radius  ar,    where 

r  =    l/(n(n-l))"1/2    .  (1.13) 


r  is  the  radius  of  the  largest  sphere  that  can  be  inscribed  into  the  unit 
simplex,  and  a  is  a  constant  such  that  0<a<l.  a  provides  a  margin 
which  ensures  that  the  algorithm  doesn't  select  a  point  outside  the 
sphere  due  to  round-off  error.  By  convexity,  an  optimal  solution  will 
occur  at  the  boundary  of  the  sphere.  Thus,  the  additional  constraint 
can  be  stated  as  an  equality.  Now  we  have, the  following  version  of  the 
LP: 
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Min 


cTD  y 


s.t.      A  D  y 

iTy 
(y-y°)T(y-y°) 


=    o 

=    1 

.     J   2 


(LP5) 


Before    the    problem    can    be    solved,     one    more    transformation 
that    moves     the     center    of     the     sphere     to     the     origin     is     useful.     Let 


-    TT     x     ,.0 


T„o   _ 


y  -  y  +  y   ,    and   eliminate    the   constant   terms   ADy°   =   0,    l1y°   =   1      and 
c    Dy    .      For  convenience,    define   B   by 


B   = 


A  D 
iT 


(1.14) 


and  the   gradient     c     of  the  objective  function  of  LP5  by 


c  =   cTD    . 


(1.15) 


Then  LP5  can  be  restated  as, 


Min 

s.t. 


-T- 

B  y 

-T— 

y  y 


=    o 

.     J  2 
=    a  r 


(LPS) 


In   order   to    solve   LP6,    note   that   an   initial   feasible    solution    to 

LP6   is      y   =    0      (which   corresponds   to   y  =   y°   in   LP5) .    This   is   also   the 

rr_  9    2 

center    of    the    sphere    y    y  =     a  v   .    An    optimal    solution    to    LP6    can    be 

obtained  by  finding  a  direction  of  maximum  rate  of  ascent  c  that  is 
feasible  with  respect  to  By  =  0,  and  moving  in  direction  -c  (maximum 
rate  of  descent)  a  distance  ar  from  y  =  0  to  the  boundary  of  the 
sphere. 

A  feasible  direction  of  maximum  ascent  is  found  by  orthogo- 
nally projecting  c  onto  the  null  space  of  B  (see  [Ref.  1:  page  17]  ), 
i.e.,    the  following  problem  has   to  be   solved   in   terms  of  c 
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Min     cTc  -   2cTc„  +  c„Tc„         =     Min      ( ||  c  -   c_J|  9)2  (1.16) 

CP  CP 

s.t.      Be  =     0   . 

Define  the  Lagrangian     L(c    ,  A.)      by 

L(cp,A)   =  cTc  -  2cTcp  *  cpTcp  -    XTBcp.  (1.17) 

First-order  optimality  conditions  for  the  Lagrangian  yield 

-2c  +  2cD  -   BTA    =  0,  (1.18) 

and 


Bcp  =  0.  (1.19) 

After  multiplying  by     B     and  dropping     2Bc    =0    (equation  1.19)   we   get 

-2Bc     =     BBT\     .  (1.20) 

T    - 1 
Assuming      (BB    )        exists,    we  can  solve  for    A 

A  =   -2(BBT)_1Bc.  (1.21) 


Substituting  equation   1.21   into  equation   1.18   gives 

c     =  c  -   BT(BBT)"1Bc.  (1.22) 

The   direction  of  maximum  rate  of  ascent  is   then 


c  =  cp/||cp||  .  (1.23) 

Thus,    the  optimal  solutions   to  LP6  and  LP5  -are 

y  =   -  *rc,  (1.24) 
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and 


y1  =  y°  -<*rc  (1.25) 


respectively. 

Finding  cp  is  the  key  part  of  Karmarkar's  method  because  it 
involves  the  major  portion  of  the  computational  work.  Solving  equation 
1.22  for  cD  can  be  viewed  as  solving  a  linear  least- squares  problem 
(see   Chapter  II.  C). 

With  the  optimal  solution  to  LP5  in  y-  space,  the  method 
proceeds  by  transforming  that  solution  back  into  x- space  by  use  of 
equation  1.12.  Next  it  defines  a  new  matrix  D  and  iterates  until  it 
reaches  the   stopping  criterion    [Ref.    1:    p.    14] 

(cTxk)/(cTx°)    <  2_c*  (1.26) 

T 

where      q      is   a   termination   parameter.    Note   that    ex   =    0   at   optimality. 

Table   1   shows  an  outline  of  Karmarkar's  basic  method. 
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Input 


Begin 


While 


TABLE  1 
ALGORITHM  1 


Problem   Size 


n 


Coefficient  Matrix. A 


Cost  Function c 

Initial  Feasible  Solution. . . x' 

Termination  Parameter q 

Feasibility  Parameter tf 


k   = 


x   = 


x 


f°  = 

cTx° 

f    = 

oo 

r      = 

l/(n(n 

-d)-1/2 

(f/f° 

>    2 

■q) 

D 

= 

diag    (x) 

y 

= 

(D"1x)/(lTD"1x)    = 

=    1 

c 

= 

cTD 

B 

= 

'    A   D    " 
.    1T       . 

CP 

= 

c    -    BT(BBT)'1Bc 

c 

= 

CP/||CP" 

y 

= 

y   -  a  re 

X 

= 

(Dy)/(lTDy) 

f 

= 

T 
ex 

k 

— 

) 

c    +    1 

End (While) 
End 
Output  Solution x 

Objective  Function  Value... f 
Iteration  Count : k 
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II.    A  VARIANT  OF  THE  PROJECTIVE  METHOD 

A.  INTRODUCTION 

The  following  sections  outline  a  variant  of  the  projective  method 
that  was  proposed  by  Wood  [Ref .  6]  .  It  has  a  practical  method  of 
bringing  an  LP  into  the  canonical  form,  and  uses  the  non-linear  objec- 
tive function  of  LP3  to  find  a  gradient  c'  which  corresponds  to  c  of 
equation  1.16.  Similar  to  the  simplex  method,  the  proposed  algorithm 
uses  a  ratio  test  to  determine  a  feasible   step  length. 

Other  practical  features  of  the  proposed  method  include  the  relax- 
ation of  the  requirement  to  know  the  optimal  objective  value  in  advance, 
and  exploitation  of  sparsity  of  A.  The  implementation  is  especially 
concerned  with  controlling  fill-in  during  the  solution  of  the  normal 
equations. 

B.  DERIVATION  OF  THE  MODIFIED  ALGORITHM 

Given  an  LP  problem  of  the  form  LP1,  we  use  a  single  artificial 
variable  x    +-.    to  attain  initial  feasibility: 


(LP7) 


Min 

(cT,M)(x,xn  +  1) 

s.t. 

Ax-(Ax°-b)xn  +  1 

= 

b 

x'xn+l 

> 

0 

where  M  is   the   cost  for  the  artificial  variable  and   x°   is   the   initial   solu- 
tion with 

(x°'x°n+l)    =   lT    •  (2-1} 

The  transformation  has  added  one  more   column  to  the  problem. 

To    get    a    homogeneous    right-hand    side    in    LP7    we    introduce    an 
additional  variable,    and  apply  the  following  projective  transformation, 
whose  inverse  is   given  by 
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(x,xn+1)   =   (n+2)(x',x'n  +  1)/x'n+2,  (2.2) 

(x',x'n  +  1)   =    (n  +  2)(x,xn  +  1)/(lT(x,xn+1)+l)  (2.3) 


x'n  +  2  =    (n  +  2)/(lT(x,xn  +  1)+l)    .  (2.4) 

Let  the  artificial  column  be  denoted  by  a,  i.e.  a  =  -(Ax°-b).  The 
transformed  problem  can  now  be  restated,  with  the  exception  of  the 
objective  function,    in  Karmarkar's  canonical  form, 


Min 


(cT,M,0)(x',x'n  +  1,x'n  +  2) 


x'n  +  2 


s.t.  (A,a,-b)(x',x,n  +  1,x,n  +  2)      =     0  (LPS) 

lT(x,'x'n  +  l'x'n  +  2)  =     n  +  2 

(x''x'n  +  l'x'n  +  2)  *     °- 


The  above  transformation  adds  yet  another  column  to  the  problem,  but 
has  the  advantage  that  it  doesn't  change  the  sparsity  of  A,  as  opposed 
to  Karmarkar's  proposal.  Rather  than  projecting  the  LP  problem  onto  a 
unit  simplex,  it  is  projected  onto  an  (n  +  2) -simplex.  This  is  done  to 
improve  numerical  stability. 

The    projective    transformation,     equation     1.12,     is    applied    to    LP8 
which   gives 


Min       (cT,M,0)D(y,yn  +  1,yn  +  2)/dn  +  2yn  +  2 


s.t.       (A,a,-b)D(y,yn  +  1,yn  +  2)      =     0  (LP9) 

iT(y,yn+1,yn+2)  =    n+2 

(y,yn+1,yn+2)  *    o  . 


The  gradient  (compare  with  equation  1.16)  of  the  objective  function  of 
LP9,  evaluated  at  the  initial  feasible  starting  point  y  =  1  is  propor- 
tional to 
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c'  =    (c1d1,...,cndnJMdn  +  1,-cTx).  (2.5) 

The    step    of    normalizing    c      in    Algorithm    1    is    replaced    by    a    ratio 
test.    Let 


jmax  =  argmax{(Cp)j>.  (2.6) 

This  allows  the  algorithm  to  make  a  step  outside  the  inscribed  sphere, 
but  maintains  feasibility  by  restricting  a  solution  to  lie  inside  the 
simplex.      Then,    the   update  of  y  becomes 

y  =  i  -  *y(cpW  >  (2-7) 

where  p  is  a  parameter  to  maintain  feasibility  such  that  0<p<l.  Table  2 
shows   the  modified  algorithm. 

C.      THE  LINEAR  LEAST- SQUARES  PROBLEM 

Computation  of  the  projected  gradient  c  during  every  iteration 
accounts  for  most  of  the  computational  workload  in  any  algorithm  based 
on  Karmarkar's  method.  Solving  equation  1.22  can  be  viewed  as  solving 
the  following  linear  least-squares  problem, 

.    Min      (||c'    -   BTAII2)2  (2.8) 

T 
where   B      is  an    (n+2)xm  matrix  with     m<(n+2)      assumed. 

T 

If   rank(B    )=m,    then   the   solution   to    (2.8)    is    given   by   the    solution 

to  the   system  of  normal  equations 

BBTA    =   Be'    .  (2.9) 

The  projected  gradient  c  is  then  the  residual  vector  of  the  least- 
squares  problem   (2.8),    i.e., 

Cp  =  C«   -   BTA     .  (2.10) 
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TABLE  2 

ALGORITHM  2 

Input 

Probl 

em  Size . . . 

n 

Coe'ff 

icient 

Mat 

rix A 

Cost 

Function. . 

0T 

Initial  Feasibl 

e  Solution. . . x° 

Termination 

Parameter t 

Feasi 

bility 

Parameter P 

Begin 

k   = 
x   = 
f  °  = 

1 
x° 

oo 

While 

(f°  " 

f(x) 

>  t 

) 

D 

= 

diag  (x) 

f° 

= 

f(x) 

y 

= 

(n+2) (D_1x)/(1TD" 

^x)  • 

c' 

= 

V(y) 

B 

= 

"AD" 
.  1T   J 

CP 

= 

c' -  BT(BBT) _1Bc' 

^max 

= 

argmax(cD)  • 

ir     J 

y 

= 

1  -  ?Cp/(Cp)imax 
(n+2)(Dy)/(iTDy) 

x 

= 

k 

= 

k  +  1 

End ( Wh 

ile) 

End 

Output 

Solut 

ion. . . 

X 

Objec 

tive  Funct 

ion  Value .  .  .  f 

Itera 

tion  Count 

k 
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T  T 

BB       is    symmetric   and   positive    definite,    given    that    Bx    has    full   column 

rank. 

As   Heath    [Ref.    7:    p.    499]    points   out,    the   ideal   choice   for   solving 

T 

the   normal   equations,    given   full   rank,    is    Cholesky    factorization.    If   B 

is  not  of  full  column  rank   the   cross-product  matrix  will  be   singular  and 

T    - 1  T  T 

(BB    )         ceases    to    exist.     BB       could    become    nearly    singular    if    B       is 

near    rank    degenerate.     Shanno     [Ref.    5:     p.     25]     shows    that    this    will 

happen    if    the    optimal    solution    is    degenerate,     i.e.     as    the    optimum    is 

approached  numerical  problems  arise   and   Cholesky   factorization  is   likely 

to  fail. 

One   problem  that   cannot  be  avoided  when   solving   the   normal  equa- 

T 
tions    is    the    fact    that    the    P-condition    number    of    BB       is    the    square    of 

that  of   BT    [Ref.    8:    p.    223],    so  that  when   BT  is   already   ill-conditioned 

it  may  be  impossible  to  find  an  accurate   solution  to  equation  2.9. 

T 
Another    important    consideration    when    computing    BB       is    that    the 

T  T 

sparsity    in    B       will    not    automatically    guarantee    sparsity    in    BB     .       In 

fact,    the    addition   of   variables    in   LP7    and    LP8   has    added    two   possibly 

T  T 

dense     bottom    rows     into     B     .       Thus,     BB        will    be     completely     dense. 

However,    one   can   cope   with   that   by   initially   omitting  the   dense   rows   in 

T  T 

B       from    the    computation    of    BB     ,     and    later    updating    the    solution    to 

equation     2.9      using     procedures        similar     to     the     ones     described     in 

[Ref.    9:    p.    58-65]. 

D.       IMPLEMENTATION 

Algorithm  2  has  been  implemented  in  FORTRAN  H  (Extended) 
Opt(2)  on  an  IBM  3033  AP  under  VM/CMS.  All  floating  point  arithmetic 
is  performed  in  double  precision.  The  program  is  designed  to  accept 
different  solution  modules  from  available  software  packages  for  solving 
the  linear  least- squares  problem. 

Input  data  sets  are  in  standard  MPS  format.  The  numerical  values 
of  the  non-zero  elements  of  the  constraint  matrix  A  are  stored  column- 
wise in  a  real  array.    For  versatility  a  full  set  of  pointers  are  defined: 

1.  IC     column  index 

2 .  R      row  index 
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3.  AP     location  of  first  non-zero  element  in  a  column 

4.  RP     location  of  last  non-zero  element  in  a  row 

5.  LINK   location  of  next  non-zero  element  (backwards) 

in  a  row. 

Brameller,    Allan  and  Hamam    [Ref .    10:    p.    104-110]    give  a   comprehensive 
discussion  of  sparse   storage   schemes. 

The  transformation  of  an  LP  problem  into  the  canonical  form  adds 
one  dense  row  and  two  dense  columns  to  the  A  matrix.  The  dense  row 
originates  from  the  convexity  constraint  equation  1.1,  the  first  dense 
column  stems  from  the  single  artificial  variable  that  is  needed  to  attain 
initial  feasibility,  and  the  second  dense  column  is  added  to  make  the 
system  of  equations  homogeneous  (see  LP8) .  The  artificial  column  is 
updated    with    the    residual   of    the    current    solution    as    long   as    the    total 

infeasibility  is  above  a   specified  threshold. 

T  T 

As    mentioned    earlier,     dense     rows     in     B        yield     BB        completely 

dense.      With   the    following  method,    this   problem   can    be   alleviated.    The 

T 
method   applies    to   any   number   of   dense    rows    in   B    ,    but    in   this    study 

we    are    only    concerned    about    two    dense    rows,     namely    the    ones    that 

result   from   the    transformations   that   are    performed   to    get   from   LP2,    to 

LP8    via   LP7.       Consider    the    projection    problem   of    equation    1.16    in    the 

following  form, 

Min      (  ||c'    -   cp  ||2)2  (2.11) 

s.t.      Be     =   0 

IT 

Replace  c      by  z,    and  let 

B     =      (B1,B2)  (2.12) 

c'      =      (c\,c\)  (2.13) 

z     -      (Zi,z2)T  -  <2-14) 
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where   B-^   is  an    (mxn-^) ,    B2   an    (mxn2)    matrix,    and  n,   is   the   number  of 
dense  columns  in  B.    Equation  2.11  now  becomes 

Min     (  ||c*2   -   z2||  2)2   +    (||  c\   -   z±\\  2)2  (2.15) 

s.t.      B1Z1      =      "B2z2 

and  when  separated,      2.15  becomes 

f   (||  c'2   -   z2||  2)2   +  Min     (||  c\   -   zJI  2)2  (2.16) 

Min     I  z-i 

2      I  s.t.    B-j^     =   -B2z2    . 

Solving  the  inner  minimization   first,    the  Kuhn-Tucker  conditions   yield 

c'l   "    zl      =   B1T^  >  (2-1?) 

and 

B1z1     =    -B2z2    .  (2.18) 

Multiplying  2.17  with  Bi    we   get, 

B1c'1   -   B1z1  =  BjB^X  .  (2.19) 

Substituting   2.18  into   2.19   gives 

B1c,1   +  B2z2  =   B1B1TA.  (2.20) 

Let    An  solve 

Bic'i  =  BiBiTA  (2-21) 

and  let    A;   solve 

(B2)j   =   BjB^A  ,    j   =   l,...,n2  (2.22) 
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where    (Bo):   is   the  j-th   column  of   Bo,    then   2.22   corresponds   to   solving 
the  matrix  equation   system 

B2  =   B1B1TA  (2.23) 

where     X    is    an    (mxno)    matrix.      Then,    the    solution    to     \    for   the    inner 
minimization  is, 

A  =    AQ  +    Az2    .  (2.24) 

Substituting  2.24  into  2.17   gives 

c'l   "    zl  =   B1T^0   +   B1TAZ2    •  (2.25) 

T  T 

To  simplify  notation,    let     h  =   Bi      Xq  and  H  =   Bi      \ 

Then  equation  2.25   becomes 

Z1  =  c\  -  h  -  Hz2    .  (2.26) 

Substituting  2.26  into  2.15,    the  overall  minimization  then  becomes 

Min      (||  c'2    -    z2U  2)2    +   (I!  h  +   Hz2||  2)2  (2.27) 

z2 

which   is   a    simple   least- squares    problem.      The   first-order   Kuhn-Tucker 
optimality  conditions  yield, 

c'2   -   hTH  =    (I   +  HTH)z2  (2.28) 

which   gives,    solving  for   z2, 

z2  =   (I  +  HTH)_1(c'2   -   hTH)    .  (2.29) 

With  the   solutions   to  2.29  and  2.26  we  have   the   desired  result, 

cp  =   z  =    (z1)Z2)T    .  (2.30) 
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In  practice  the  following  procedure  is  followed: 


1. 

2. 
3. 
4. 
5. 
6. 
7. 
8. 

9. 

10 
11 


compute  BjB-i 

T 
factor  BtBt  ,  e.g.  using  Cholesky  factorization 

rn  m 

solve  B-j^B-^  A  =  B^c'^,  the  solution  is  AQ 

T 
compute   h  -  Bi  A Q 

solve   BiB]_  A  =  (Boji,  the  solution  is   A-, 

m 

solve   BiBi  A  =  (B2)2/  tne  solution  is   A2 

T 

compute   H  =  B-,  }±. 
compute    ( I 


+   HTH)-1 


(note   that  HXH   is 


2x2 


matrix  if  there  are  2  dense  rows  in  B  ) 

rn 

compute   (c'a  -  h  H) 
compute  Zo 
compute  Zi 


The  given  procedure  is  efficient  since  the  factorization  of  BB^  is 
computed  only  once  and  the  same  system  is  solved  three  times  using 
different  right  hand- sides  each  time. 

As  a  stopping  criterion  for  the  algorithm  the  following  rule  is 
used, 


IF     argmax    ( |  x,k   -   x^"1!)    <  t     STOP   , 
j  J  J 


where  t  is  a  real  constant.      The  convergence  criterion 


(2.31) 


(||  cp  ||/|cTx°|)   <  t 


(2.32) 


mentioned   by   Lustig    [Ref.    4:    p.    12]    can   also   be   used   to   terminate    the 
algorithm. 


25 


III.    SOLVING  THE  LINEAR  LEAST- SQUARES   PROBLEM 

A.  INTRODUCTION 

In  this  study,  the  linear  least- squares  problem  is  solved  by 
explicitly  computing  and  solving  the  normal  equations.  Although  the 
normal  equation  approach  experiences  problems  when  applied  to  ill- 
conditioned  or  near  rank  deficient  matrices,  it  behaves  acceptably  with 
sparse  and  well- conditioned  matrices    [Ref .    11] . 

The  numerical  methods  for  solving  normal  equations  generally  fall 
into  two  classes,  direct  and  iterative  methods.  One  representative  of 
each  class  is  considered  in  this  study.  Little  can  be  said  as  to  which 
class  of  methods  is  better,  except  that  direct  methods  are  more  attrac- 
tive in  terms  of  computational  work,  and  iterative  methods  may  require 
less   storage    [Ref.    12:    p.    11]. 

B.  CHOLESKY  FACTORIZATION 

The  method   implemented  is   given   in   George   and  Liu    [Ref.    12],    and 

uses    Cholesky    factorization   with    a    minimum- degree    ordering    to    solve    a 

T 

large    sparse    positive    definite    system    of    equations.     Since    BB       is    only 

guaranteed  to  be  positive  semidefinite,  a  modification  to  the  Cholesky 
factorization  algorithm  is  considered  in  Chapter  IV.  A.  3.  to  accommodate 
the   semidefinite  case. 

The  minimum-degree  algorithm  is  a  reordering  heuristic  which 
attempts  to  reduce  fill-in  during  the  factorization  phase.    The  reordering 

phase    is    entirely    symbolic;    it    amounts    to    a    symmetric    row    and    column 

T 
permutation    of    BB       which    corresponds    to    reordering    the    columns    in 

t  nr 

B     .    During    this    phase,    BB       doesn't    have    to    be    computed    numerically; 

only    its    structure    has    to    determined.     Also,     the    factorization    is    first 

performed    symbolically,     thus    allowing    a    static    data    structure    for    the 

Cholesky   factor   L.      An   outline    of   the   phases    of   the   algorithm   is    given 

in  Table  3.    See  also  Heath    [Ref.    7:    p.    499]. 
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TABLE  3 

MINIMUM- DEGREE  ORDERING  /  CHOLESKY  FACTORIZATION 

ALGORITHM  (MDOC) 

1. 

Determine  the  nonzero  structure  of  BB  . 

2. 

Find  a  permutation  matrix  P  such  that  PBBTPT 

has  a  sparse  lower  triangular  Cholesky  factor 

L. 

T  T 
Factor  PBB  P   symbolically  and  set  up  the  data 

3. 

structure  for  L. 

4. 

T  T 

Compute  PBB  P   numerically. 

5. 

Factor   PBBTPT  =  LLT  numerically. 

6. 

Solve   Lz  =  PBb '   (back  substitution). 

i : 

T 
Solve   L  y  =  z   (forward  substitution). 

8. 

x  =  PTy. 

1.        The  Minimum- Degree  Ordering  Heuristic 

The    heuristic    finds    an    ordering    of    a    symmetric    matrix    such 

that   fill-in   is   low  when   the   matrix   is   being   factored.    The   basic   idea   is, 

T 
at    each     (simulated)     factorization    step,     to    permute    the    part    of    BB 

remaining  to  be  factored   so  that  a  column  with  the  fewest  nonzeros  is   in 

the    pivot    position.    The    implementation    consists    of    six    subroutines    that 

are    given   in    George   and   Liu    [Ref.    12:    pp.    124-137].      The    subroutines 

accept    as    input    the    adjacency    graph    associated    with    BB       represented 

by   an   adjacency   structure,    and   return   as   output  a    symmetric   permuta- 

T  T 

tion  of  BB      given  as  a  permutation  vector  for  the  columns  of  B     . 

T 
Let      G=(X,E)       be    the    adjacency    graph    of    BB     ,     where    X    is 

the  set  of  nodes,  and  E  is  the  set  of  edges.  Then,  the  nodes  corre- 
spond to  the  variables  of  the  least- squares  problem,  i.e.  the  columns  of 
B  .  Two  nodes  x  and  y  are  said  to  be  adjacent  if  {x,y}  is  an  edge  in 
E.    The   adjacent   set  of  Y,    Y<=X  is   defined  and   denoted  by 
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Adj(Y)    =    (x    e    X-Y|{x,y}eE     for  some  yeY).  (3.1) 

An  adjacency  list  for  X£X  is  a  list  of  all  nodes  in  Adj({x}). 
Finally,  an  adjacency  structure  for  the  graph  G  is  the  set  of  adjacency 
lists  for  all  xeX  [Ref.  12:  pp.  37-41].  The  particular  adjacency 
structure  used  in  the  ordering  heuristic  stores  elements  in  each  adja- 
cency list  in  contiguous  locations.  An  entry  point  array  to  the  first 
element  in  each  list  allows  access  to  the  list. 

The  ordering  heuristic  is  based  on  graph  theory,  and  involves 
the  notion  of  elimination  graphs,  quotient  graphs,  reachable  sets  and 
indistinguishable  nodes.  George  and  Liu  [Ref.  12:  pp.  92-124]  may  be 
consulted  for  more   details. 

2.       Factorization  and  Solution 

The  components  of  the  lower  triangular  Cholesky  factor  L  of 
BB  are  computed  using  the  so  called  "inner  product  form"  algorithm 
[Ref.    12:    p.    20].      The  elements  of  L  are   given  by, 

'ji=    (ejj   ■%    l2Jk>V2         for  j=l,2,...,m  .  (3.2) 


ljj  =    (ey   -   2     Wjk'^jj         for  i=J+l.J+2 m  (3.3) 

where  the     e-   are  the  elements  of  BB     . 

After    the    factorization    has    been    computed,     the    following    two 
linear  systems  have  to  be   solved    (see  also  Table  3) : 


Lz  =   Pb'    ,  (3.4) 


and 


LTy  =   z   .  (3.5) 
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Solving    system    3.4    by    back    substitution    involves    the    use    of    "inner 
products"    [Ref.    12:    p.    25],    defined  by 

i-l 
Zj  =    (b'j   -   2     lik^/lii  for  i=l,2,...,m     ,  (3.6) 

K  -  1 

where    b'j    stands   for   the    i-th   element   of    the    right-hand    side    of    (3.4). 
For  this  case  L  must  be  accessed  row-by-row. 

System    3.5    is     solved    by    forward    substitution    using     "outer 
products"    [Ref.    12:   p.    26],   defined  by 

yj  =   Zj/ly         for  i=l,2,...,m  (3.7) 

(zi+1,..,zm)  <-  (zi+1)..,zm)  -  yiCVi^,..,^!)  . 

T 
For  the   latter   case,    L      must   be   accessed   r-ow-by-row,    or   L  column-by- 
column  instead. 

C.       INCOMPLETE  CHOLESKY  FACTORIZATION 

This  method  is  an  implementation  by  Ajiz  and  Jennings  [Ref.  13]  of 
the  incomplete  Cholesky  conjugate  gradient  algorithm  (ICCG),  whose 
theory  is  given  in  [Refs.  14,15].  The  algorithm  requires  that  the 
coefficient  matrix  of  a  set  of  simultaneous  linear  equations  be  symmetric 
and  positive  definite.  It  consists  of  two  distinct  parts,  one  being  the 
Cholesky  factorization,  which  can  be  complete  or  incomplete,  and  the 
other  a  conjugate  gradient  iteration  to  solve  a  preconditioned  linear 
system,    where   the   Cholesky  factor   serves  as   the  preconditioner. 

Golub  and  Van  Loan  [Ref.  16:  pp.  373-377]  point  out  that  precon- 
ditioning is  essential  for  obtaining  good  convergence  rates  with  conju- 
gate   gradient    methods.     The    convergence    rate    is    closely    linked    to    the 

P-condition    number,     which    is    the    ratio    of    the    largest    to    the    smallest 

T 
eigenvalue.    It    was    mentioned    earlier    that    the    condition    number   of    BB 

T 
will    be    the    square    of    the    one    of    B     .       Ill-conditioned    problems    have 

large    condition    numbers,     and    hence    slow    convergence.     Preconditioning 

is    a    process    of    transforming    a    linear    system    so    that    its    P-condition 

number  is  improved    [Ref.    17:    p.    979]. 
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Consider  again   the   original  linear  least-squares   problem    (2.9)    with 
a  generic  right-hand   side  b'  and  x  the  unknowns, 

BBTx  =  b'.  (3.8) 

Then,    equation    3.8    can    be    preconditioned    with   a    transformation   matrix 
L  giving 

L_1BBTL"Tx  =  I/V,  (3.9) 


or 


L"1BBTL"Ty  =  b,  (3.10) 


where     y     =     L    x     and     b     =     L     b'.        According     to     Ajiz     and    Jennings 

[Ref.    13:     p.     950]     the    ideal    choice    of    transformation    matrix    L    is    the 

T 
Cholesky  factor  of  BB     ,    since 

L_1BBTL"T  =   I    ,  (3.11) 


provided  one   could  perform  exact  arithmetic. 

The    objectives    of    the    incomplete    factorization    phase    of    the    algo- 

T 
rithm    are    to    transform    BB       as    close    as    possible    to    I,    and    to    reduce 

fill-in  in  the  factor  L.  This  is  accomplished  by  discarding  some  off- 
diagonal  coefficients  during  the  factorization,  whose  magnitudes  fall 
below  a  preset  threshold  limit.  The  result  of  this  operation  is  an  incom- 
plete  Cholesky  factor  L  that  must   satisfy 

BBT  =  LLT   -   C    ,  (3.12) 

where  C  is  the  matrix  of  elements  omitted  from  the  factorization . 
Unfortunately,  omission  of  elements  from  the  factorization  process  can 
destroy  the  positive  definiteness  property  and  hence  lead  to  a  break- 
down   of    the    process.    Ajiz    and    Jennings    [Ref.    13:    pp.    950-951]    have 
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shown  that  introducing  diagonal  modifications  in  C  will  retain  the  posi- 
tive definiteness  property.  The  matrix  C  will  be  symmetric  and  will 
have   diagonal   elements   that  are   greater  than  or  equal  to  zero.    Thus,    C 

is   a   positive   semidefinite   matrix.    Assuming   BB      to   be   positive   definite, 

T 
adding  C  will  result  in  LL      also  being  positive   definite. 

Convergence  of  conjugate  gradient  iterations  is  only  guaranteed 
for  the  positive  definite  case.  Thus,  a  modification  to  adapt  the 
Cholesky  factorization  to  the  positive  semidefinite  case  as  with  the 
direct  method  may  cause  slow  convergence.  The  modification  is  consid- 
ered in  Section  IV.  A.  3.  Table  4  gives  an  outline  of  the  ICCG 
algorithm. 


TABLE  4 

INCOMPLETE  CHOLESKY       CONJUGATE  GRADIENT   (ICCG) 

ALGORITHM 

1. 

Obtain   L,    an   incomplete   Cholesky 

factor    of    BBT. 

2. 

Solve         Lb   =   b'     for   b   by   forward 

substitution. 

3. 

Solve         L"1BBTL"Ty      =      b      for      y 
gradient    iteration. 

by      conjugate 

4. 

Determine   x   by  back    substitution 

T 

in      L  x   =   y. 

1.       The  Incomplete  Factorization 

The    procedure    presented    here    is    given    in    Ajiz    and  Jennings 

[Ref.    13:      pp.      951-952],      and      Jennings      and      Malik      [Ref.  14:      pp. 

310-313]  .       To    see    how    the    elements    in    column    j     of    L    are  computed 

consider    the    following.     From    matrix    equation    3.12    a    typical  elemental 
equation  may  be  written  as 

Vij   =   Sij    *   °ii   "  £    Wjk  <*j)  (3.13) 
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where   the   subscripts   for  elements  1  refer  to  their  positions   in  the  lower 

T 
triangular  factor  L,    and  the  e-  are  the  elements  of  BB     .      Let 


3-1 


e 


»  =  e« " &  Vik  •  -  <»■"> 


Then,  assuming  the  elements  of  L  have  been  computed  for  columns  1  to 
(j-1),  all  the  e  ■■  in  column  j  can  be  computed.  First  consider  the  case 
where  all  the  e'--  pass  the  rejection  test,  i.e.,  are  to  be  retained. 
Hence,    by  setting  c-pO  in  equation  3.13,    we   get 

'ij  •  -Vji  •  (315) 

In  case  any  e' •■  is  to  be  rejected,  1-  is  set  to  zero,  implying  c-=-e'-. 
The  off-diagonal  term  c-  implies  diagonal  additions  c-  and  c-  to  the 
matrix  C  in  order  to  have  the  positive  semidefinitness  property.  The 
matrix  C  doesn't  have  to  be  stored  in  memory;  only  the  diagonal  addi- 
tions c--  are  of  further  interest. 

Before  any  1-  can  be  computed,  L-  has  to  be  determined.  With 
all  the  diagonal- additions  c--  that  resulted  from  rejections  of  e  '•■  -during 
the  computation  of  columns   1   to   (j-1),    an  expression  for  L-   becomes 

I,,  =    (e*«   ♦     5     c./k))1/2  (3.16) 

where  e '  ••   is   defined  by 

eVeiJ+£  cij(k)  -  %  'ik'  '  <317> 

and  c^'  '  is  the  diagonal  addition  to  c-  resulting  from  deletion  of  eV-, 
and  c--^    '    the  diagonal  addition  to  c-  resulting  from  deletion  of  e   -^. 

The  rejection  operation  tests  the  magnitude  of  an  element  e '  •• 
in  relation  to  the  current  values  of  the  corresponding  diagonal  elements 
e-  and  e-   respectively,    whose  values  are   given  by 

"ii  =  '*jj  \?,0Jj(k)  -  (318) 

ic  *  ]+l 
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where    the    index   k   refers    to   the   rows    for   which   rejections    in    column   j 

have  had  a  bearing  on  e-     and 

JJ  j 


ii  =  eii  +A    cii(k)    •  (3-19) 


In  equation  3.19,  k  refers  to  the  columns  for  which  rejections  in  row  i 
have  had  a  bearing  on  e^.      An  element  e'--   is   rejected  if 

eY <  *2  ^ii  •  <3-2°) 

where  \p  is  the  preset  rejection  parameter.  A  choice  of  \p-0  will  retain 
all  elements,  thus  leading  to  a  complete  Cholesky  factorization.  A  choice 
of  \p-l  will  cause  all  off -diagonal  elements  to  be  rejected.  Ajiz  and 
Jennings  [Ref.  13:  p.  952]  recommend  that  the  rejection  parameter  be 
in  the  range  0.01<i^<0.2  for  effective  incomplete  factorizations. 

The  diagonal  modifications  in  C  which  result  from  the  rejection 
of  an  element  e   -  are   given  by 

cJi(k)  =  lAjKVekk)1'2  -  •      <3-21> 

and  by 

cU(k)    =    l<Mkl(eU/ekk)1/2    ■  (3.22) 

With  the  successive  application  of  equations  3.14  in  conjunction  with  the 
rejection  operation,  3.17,  3.16  and  3.15,  all  elements  in  column  j  of  L 
are  determined. 

2.       The  Conjugate  Gradient  Iteration 

The  conjugate  gradient  method  of  Hestenes  and  Stiefel 
[Ref.  18]  is  applied  to  matrix  equation  3.10.  It  uses  the  following 
vectors,    the  letter  k  indicates   the  k-th  iteration, 

a)  p(*)  conjugate  gradient  vector 

b)  r'  '  residual  vector 
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c)  y(  '  solution  vector 

d)  u(k)  product  of  L_1BBTL"T  and  p<k)  . 

The  initial  values  for  k=0  are  y(°)    =   0  and  p(°>    =   r^°)    =  b.      The  algo- 
rithm for  one  iteration  is  as  follows 

u(k)  =    L-lBBTL"Tp(k) 

ak  =   (pWjVk^pClOjT^k) 

y(k+l)  =  y(k)   +  akP(k) 

r(k+l)   =  r(k)    .  <*u(k) 

(3  _.   (r(k+l))Tr(k+l)/(r(k))Tr(k) 

p(k+l)   =  r(k+l)    +  ^(k)    i 

The   first   step   in   the   above   algorithm   is   obtained  without   computing   the 
transformed  matrix  explicitly  by  the  following  three  operations, 

1.  LTv(k)    =   p(k)  (back    substitution) 

2.  w(k)    =    BBTv(k)       (pre-multiplication) 

3.  Lu(k)    =   w^k)  (forward    substitution). 

The  algorithm  is  terminated  when 

llr(k)||  /  ||b||  <  tolerance  (3.23) 

where   b  equals   the    starting   residual   r^   ',    since   y^    '    was    chosen   to  be 
zero. 
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IV.    MODIFICATIONS  AND  COMPUTATIONAL  RESULTS 
A.    "ALGORITHMIC  MODIFICATIONS 

1.  Iterative  Improvements 

A  provision  to  improve  the  solutions  to  equations  2.21  and 
2.22  has  been  implemented,  because  their  residuals  are  often  unduly 
high.      The  residuals  are  computed  as 

r  =  B^^A  -  b'   ,  (4.1) 

where  A    stands    for  Aq>   Ai    and     A2    respectively    and    b'    is    generic    for 

B1Tc1',    (B2)1  and   (B2)2. 

-  fi 
If     |r-|>    g     for     some     i,     where   e     is     set    to,     say,     10      '     the 

following   systems  are   solved  for    A' 

B^^A'   =    -r.  (4.2) 

An  improved  solution  is  then  obtained  by  adding    A  and    A' 

B1B1T(A+  A')  =  r+b'+(-r)  =  b'.  (4.3) 

The  improvement  can  be  repeated  if  the  residuals  of  equation  4.3  are 
still  found  to  be  too  high.  'With  this  modification  the  direct  method 
(MDOC)    has   become  a   semi-iterative  method. 

2.  Removal  of  the  Artificial  Column 

The  update  of  the  artificial  column  with  the  residuals  of  the 
current   solution 

A(x',x'n  +  1)    -  bx'n  +  2  =  r  (4.4) 

Ax<k)    -  b  =  Ax(k_1)    -  b  -   r/x'n  +  2(k)  (4.5) 
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on  each  iteration  has  been  augmented,  such  that  the  artificial  variable 
can  be  deleted  from  the  problem  when  the  infeasibility  becomes  small. 
After  the  deletion,  only  one  dense  column  remains  in  B,  and  the  solu- 
tion of  the  least-squares  problem  is  simplified.  The  motivation  for  this 
modification  is  to  avoid  some  of  the  numerical  problems  that  arise  from 
variables  that  approach  zero. 

...  T 

3.  Positive  Semidefinite  BB 

It    can    be    shown   that    if,    and   only    if   a    diagonal   element   ever 

T 
goes    to    zero   in    a    Cholesky    factorization,    BB       is    positive    semidefinite, 

not  positive  definite.  Furthermore,  all  elements  below  the  zero  diagonal 
element  must  also  be  zero  and  a  nonunique  solution  to  the  normal  equa- 
tions can  be  obtained  by  setting  the  corresponding  A;  to  0.  This  solu- 
tion can  be  obtained  by  replacing  the  zero  diagonal  element  with  any 
positive    value    and    continuing    with    the    factorization.       Restated,     if    a 

diagonal    element    L-    in    the    partially    computed    Cholesky    factor    is    less 

-6 
than    or    equal    to    10      ,    set    1«=1,    and    1«.=0    for    i=j  +  l,  j  +  2,  .  .  .  ,m    in    equa- 
tions   3.2    and    3.3    respectively.      This    procedure    is    also    useful    to    deal 
with  numerical  problems  which  arise  from  degeneracy  near  optimality. 

4.  Weighted  Homogeneity  Variable 

The  initial  solution  (y0,y°n  +  i  ,y0n  +  o)  ~  1  used  to  begin  the 
projective  algorithm  is  arbitrary,  and  in  some  sense,  the  "homogeneity 
variable"  y°n  +  9  *s  fundamentally  different  than  the  other  variables. 
Consequently,  a  modification  has  been  made  to  allow  weighting  the 
starting  solution  y°n  +  2  differently  from  the  other  y°-,  i=l,  .  .  .  ,n  +  l.  Let 
y°n  +  2=s.    With   lTy°=n  +  2,    the  other  y0,  are   set  equal, 

y°i  =    (n+2-s)/(n+l),      i=l,...,n+l.  (4.6) 

It  is  hoped  that   such  a  weighting  might  lead  to  lower  iteration  counts. 
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B.   TEST  PROBLEMS  AND  COMPUTATIONAL  RESULTS 

At  first,  a  small  test  problem  (TEST1)  was  used  to  verify  the 
correctness  of  the  source  code.  The  two  algorithms  were  then  tested  on 
seven  problems  which  have  been  also  used  for  testing  by  Lustig 
[Ref.  4:  p.  19].  Table  5  shows  some  of  the  characteristics  of  the  test 
problems . 


TABLE  5 

TEST  PROBLEMS 

Name 

Rows 

Columns 

Logicals 

Nonzeros 

Density 

TEST1 

6 

2 

6 

10 

0.833 

AFIRO 

27 

32 

19 

95 

0.  103 

ADLITTLE 

56 

97 

41 

522 

0.096 

SHARE2B 

96 

79 

83 

901 

0.  119 

ISRAEL 

174 

142 

174 

2529 

0.  102 

BRANDY 

193 

249 

54 

2204 

0.046 

E226 

223 

282 

190 

2578 

0.041 

BANDM 

305 

472 

0 

2659 

0.018 

Table  6  summarizes  the  computational  results.  All  CPU  times 
represent  the  time  in  seconds  to  set  up  the  major  part  of  -the  data 
structure,  solve  the  LP,  and  write  out  a  few  parameters  on  each  itera- 
tion and  the  solution.  Times  to  read  in  the  data  from  MPS  format  are 
not  included. 

The  convergence  criterion  equation  2.31  is  used  with  t=0.05  for  all 
test  problems  but  TEST1  and  SHARE2B,  where  t=0.001  and  t=0.01 
respectively.  No  solutions  have  been  obtained  for  the  problems  AFIRO, 
BRANDY  and  BANDM:      the  algorithm  will  not  converge  to  the  optimum. 

The  feasibility  parameter  p  (equation  2.7)  is  set  to  0.9995,  giving 
the   best    overall   performance    of   the    algorithm.    Tests    with  p=0. 9999    and 
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p=0. 99999  on  the  data  sets  TEST1  and  SHARE2B,  indicate  that  with 
these  higher  values  the  number  of  iterations  necessary  to  attain  feasi- 
bility is  reduced,  but  the  total  number  of  iterations  remains  about  the 
same.    In  addition,    the  accuracy  of  the   solutions  deteriorates. 

An  intermediate  solution  is  declared  feasible,  i.e.  the  artificial 
variable  is  removed  from  the  problem,  as  soon  as  x  +iM  is  less  than 
10"  .  The  choice  of  10  is  arbitrary,  and  depends  on  the  value  of  the 
optimal  solution.  Using  alternate  values  of  1.0  or  10  makes  little 
difference  in  the  number  of  iterations  necessary  to  attain  feasibility. 


TABLE  6 

COMPUTATIONAL 

RESULTS 

MDOC 

ICCG 

Problem 

Iterations 
feasible  optimal 

Time 

Iterations 
feasible  optimal 

Time 

TEST1 

1                 6 

0.03 

1                 6 

0.04 

ADLITTLE 

13              21 

1.96 

13              21     . 

2.18 

SHARE2B 

6              27 

5.04 

7              37 

5.31 

ISRAEL 

19              98 

311 

19              89 

347 

E226 

8              49 

107 

9               49 

130 

Factorization    failures    plagued    both    algorithms    before    the    modifica- 

T 
tion    to    accomodate    a    semidefinite    BB       was    implemented.    These    failures 

occured  near  the  optimum,  e.g.  E226,  or  with  the  ICCG  algorithm  when 
setting  the  rejection  parameter  to  a  value  greater  than  zero.  After 
implementing  the  semidefinite  modification  these  factorization  problems 
have  been  cured,  but  the  ICCG  algorithm  now  shows  very  slow  conver- 
gence. For  example,  a  solution  to  SHARE2B  is  obtained  only  after 
112.36  CPU  seconds  with  ^=0.015  and  all  other  parameters  unchanged. 
Thus,    because   storage   is   not  at  a   premium,    only  complete   factorizations 
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are  used.  Satisfactory  numerical  results  with  the  incomplete  factoriza- 
tion are  obtained  only  with  the  trivial  test  problem  TEST1;  computation 
times  are  always  inferior. 

The      minimum -degree      ordering      works      very      well      in      practice. 
Computational     cost     seems     to     be     moderate,      e.g.,      0.25     seconds     for 

SHARE2B,     and    3.12    seconds    for    BANDM.     Table    7    gives    densities    of 

T 
BB      and  the   Cholesky  factors  with  and  without  reordering.      Substantial 

reductions     in     storage     requirements     were     achieved     when     doing     the 


reordering 


TABLE  7 

DENSITIES 

Re-ordered 

Not  re- 

■ordered 

Nonzeros 

Nonzeros 

Nonzen 

3S 

Name 

in   BBT 

Density 

in  L       Density 

in   L 

Density 

TEST1 

20 

0.952 

20 

0.952 

20 

0.952 

AFIRO 

62 

0.177 

107 

0.305 

194 

0.553 

ADLITTLE 

384 

0.241 

411 

0.258 

816 

0.512 

SHARE2B 

871 

0.187 

1021 

0.219 

1134 

0.243 

ISRAEL 

11227 

0.737 

11433 

0.751 

13743 

0.903 

BRANDY 

2853 

0.  152 

3429 

0.183 

9760 

0.521 

E226 

2823 

0.113 

3639 

0.  146 

10735 

0.430 

BANDM 

3724 

0.080 

4660 

0.  100 

32090 

0.688 

The    densities    and   number 

of   nonzeros    in   L   or 

T 
BB       are   relative 

to    a    symmetric    half 

of    a    matrix;    the 

diagonal 

elemenl 

s    are    not 

included. 

Computational  results  indicate  that  iterative  improvements  are  not 
always  an  absolute  necessity.  Residuals  of  a  current  solution  tend  to  be 
high    at    the    start    of    iterations,     probably    due    to    a    less    than    optimal 
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choice  of  M.  Doing  a  few  (2  or  3)  iterative  improvements  at  this  stage 
stabilizes  the  computations  until  the  algorithm  gets  close  to  the  optimal 
solution.  Numerical  problems  near  the  optimum  are  so  severe  that  even 
allowing  a  prohibitively  high  number  like  25  iterative  improvements  has 
no  apparent  influence  on  the  quality  of  the   solution. 

The  mildest  form  of  numerical  difficulty  near  the  optimum  is  slow 
convergence,  e.g.  ISRAEL  and  E226.  For  SHARE2B,  a  high  quality 
solution  is  obtained  with  no  numerical  difficulty.  This  data  set  is  chosen 
to  test  the  weighted  homogeneity  variable  modification.  S  is  given 
several  values  in  the  range  1  to  100.  Iteration  counts  range  from  25  to 
37.  The  quality  of  the  solutions  is  sometimes  reduced,  however.  Only 
s=35  (27  iterations)  and  s=50  (28  iterations)  give  high  quality  solutions 
with  low  iteration  counts.  Values  of  s  greater  than  75  creat  conver- 
gence failures. 

C.       CONCLUSIONS 

The  low  iteration  counts  for  some  of  the  test  problems  are  prom- 
ising, although  CPU  times  seem  to  tell  the  difference.  These  high  CPU 
times  result  from  test  problems  having  slow  convergence  near  the 
optimum,  which  is  believed  to  be  due  to  many  variables  going  to  zero. 
Thus,  a  technique  to  drop  variables  going  to  zero  from  the  LP  could 
well   speed  up   convergence. 

The  ICCG  algorithm  does  not  perform  very  well  on  the  test  prob- 
lems. Computational  results  are  generally  inferior  to  those  obtained 
with  the  MDOC  algorithm.  Thus,  no  further  research  into  this  method 
seems  warranted. 

In  view  of  the  numerical  problems  encountered  when  solving  the 
least- squares  problem  with  the  normal  equations  approach  and  Cholesky 
factorization,  another  method  that  does  not  use  square  roots  should  be 
considered  for  implementation.  A  very  promising  candidate  in  this 
respect  is  Givens  rotation.  See  Gentleman  [Ref.  19]  and  George  and 
Heath   [Ref.    20]. 
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